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ABSTRACT 

We  quantify  dynamieal  ehange  in  nonlinear  time- 
serial  data  via  dissimilarity  measures  between  statistieal 
distribution  funetions.  These  new  measures  are  superior  to 
traditional  nonlinear  measures,  and  give  robust  and  timely 
forewarning  of  normal-to-abnormal  transitions  in 
physiologieal  regimes.  These  novel  measures  also  provide 
rapid  quantitative  assessment  of  physiologieal  ehange, 
sueh  as  alertness,  soldier  readiness,  and  bioagent  exposure. 

1.  INTRODUCTION 

Charaeterization  of  ehange  in  eomplex  system 
dynamies  is  partieularly  diffieult  for  physiologieal 
proeesses,  whieh  share  several  eonfounding  features:  non- 
stationarity,  nonlinearity,  multiple  time  seales,  and  strong 
sensitivity  to  environmental  perturbations.  The  proeess 
dynamies  ean  be  quantified  by  traditional  nonlinear 
measures,  sueh  as:  (i)  the  first  minimum.  Mi,  in  the 
mutual  information  funetion,  whieh  is  a  measure  of 
nonlinear  de-eorrelation  time;  (ii)  the  eorrelation 
dimension,  D,  whieh  quantifies  dynamieal  eomplexity; 
and  (iii)  the  Kolmogorov  entropy,  K,  whieh  is  an  indieator 
of  predietability  (or  equivalently,  the  rate  of  nonlinear 
information  loss).  Kolmogorov  entropy  and  eorrelation 
dimension  are  usually  defined  in  the  limit  of  zero  seale 
length.  However,  all  real  data  have  noise  and  even 
noiseless  model  data  are  limited  by  the  finite  preeision  of 
eomputer  arithmetie.  Finite-seale  values  of  K  and  D  do 
not  eapture  the  full  dynamieal  eomplexity  of  these 
proeesses  in  the  absenee  of  noise  and  have  smaller  values 
than  expeeted  for  the  zero-seale-length  limit.  In  general, 
traditional  nonlinear  measures  ean  distinguish  between 
elear-eut  regular  and  ehaotie  dynamies.  However,  they  are 
not  suffieiently  sensitive  to  diseem  between  slightly 
different  ehaotie  regimes,  espeeially  when  data  are  limited 
and/or  noisy.  This  laek  of  diserimination  arises  from 
averaging  over  the  global  dynamies,  whieh  erases  most  of 
the  dynamieal  details. 

We  address  these  limitations  in  traditional  nonlinear 
measures  via  new  measures  to  quantify  ehange  in  time- 
serial  data  (ineluding  maehines)  as  follows.  We  eonvert 
time-serial  data  to  a  diserete  geometrie  (phase  spaee) 
representation.  A  distribution  funetion  deseribes  the 
visitation  frequeney  and  sequenee  of  the  diserete  phase- 


spaee  states;  (un)ehanging  dynamies  lead  to  an 
(un)altered  distribution  funetion.  Dissimilarity  measures 
quantify  the  distanee  between  a  test  ease  distribution 
funetion  and  the  baseline  distribution  funetion.  Large 
dissimilarity  means  that  the  system  is  far  from  the 
baseline,  as  a  forewarning  of  an  abnormal  event, 
partieularly  in  human  data.  A  eomparison  of  the  results 
shows  a  signifieant  and  eonsistent  superiority  of  our  new 
measures  over  traditional  nonlinear  measures  for 
deteetion  and  forewarning  of  eondition  ehange  in  real 
physiologieal  data. 

We  organize  the  paper  as  follows.  Seetion  2  diseusses 
the  analysis  methodology  and  our  new  dissimilarity 
measures.  Seetion  3  presents  the  results,  as  applied  to 
physiologieal  time  series.  We  summarize  the  eonelusions 
in  Seetion  4,  with  aeknowledgments  in  Seet.  5. 

2.  ANALYSIS  METHODOLOGY 

We  begin  with  a  proeess-indieative  sealar  signal,  e, 
that  is  sampled  at  equal  time  intervals,  x,  starting  at  time 
to,  yielding  a  time  series  e,-  =  e{to  +  h),  i  =  1,  2,  N. 
Physiologie  data  are  usually  affeeted  by  artifaets  (e.g.,  eye 
blinks,  musele  twitehes).  We  remove  essentially  all  of 
these  artifaets  by  fitting  a  parabola  to  e,  over  a  moving 
window  of  2w-l-l  points  with  the  same  number  of  data 
points,  w,  on  either  side  of  a  eentral  point.  We  use  the 
fitted  value  at  the  eentral  point  to  estimate  the  artifaet 
signal,  f.  Subsequent  analysis  uses  only  the  artifaet- 
filtered  data,  v,  =  e,  -f ,  whieh  is  essentially  artifaet-free. 

We  extraet  the  main  features  of  the  underlying 
dynamies  via  a  dynamieal  reeonstruetion  teehnique. 
Speeifieally,  we  form  a  (7-dimensional  time -delay  veetor, 
y{i)  =  [Xi,  Xi+x,  .x:/+(d-i):t].  The  ehoiee  of  lag,  L,  and 
embedding  dimension,  d,  determine  whether  the  data  are 
under-sampled  (projeeted  upon  themselves)  or  redundant. 
This  phase-spaee  approaeh  is  used  for  both  our  new 
measures  of  eondition  ehange  and  the  traditional 
nonlinear  measures;  see  the  referenee  for  more  details. 
We  also  note  that  different  system  observables  eontain 
dissimilar  information.  Thus,  phase-spaee  reeonstruetion 
ean  be  easier  from  one  variable,  but  more  diffieult  (or 
even  impossible)  from  another.  Our  analysis  seeks  a 
balanee  between  these  eaveats  and  the  praetieal  eonstraint 
of  a  limited  amount  of  noisy  data. 
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Sensitive  change  discrimination  uses  the  phase-space 
reconstruction  of  the  reconstructed  dynamics,  as  follows. 
We  first  represent  each  continuous  value,  as  a  discrete 
symbol,  s,,  that  is  one  of  S  different  integers,  0,  1,  . . .,  5-1 : 

0  <  ^ ,  =  INT[S{Xi  -  {x„ax  -x„,i„)\<S-\.  ( 1 ) 

The  function,  INT,  converts  a  decimal  number  to  the 
closest  lower  integer.  The  minimum  and  maximum  values 
of  Xi  are  Xmi„  and  x^ax,  respectively,  over  the  baseline 
(reference  data).  We  require  that  Si(Xi  =  Xmax)  =  5  -  1  in 
order  to  maintain  exactly  S  distinct  symbols.  The 
discretized  vector,  y{i)  =  [^„  Si+i,  ...,  partitions 

the  phase  space  into  S‘‘  hypercubes  or  bins.  We  count  the 
number  of  phase-space  points  in  each  bin  to  obtain  a 
distribution  function.  We  denote  the  population  of  the  i-th 
bin  of  the  distribution  function  as  Qi  for  the  baseline,  and 
Ri  for  a  test  case.  The  parameters  (5,  N,  d,  and  depend 
on  the  specific  data,  and  reflect  the  underlying  dynamics. 

We  compare  the  distribution  functions  of  the  test  case 
and  baseline  dynamics  by  measuring  the  difference 
between  g,  with  7?, by  the  statistics  and  Li  distance: 


x^=z(e,-R,)^/(a+R^, 

i 

(2) 

L  =  m-Ri 

(3) 

The  sums  in  Eq.  (2  -  3)  run  over  all  of  the  populated 
phase-space  cells.  These  measures  account  for  changes  in 
the  geometry  and  occurrence  frequency  of  the  distribution 
function.  The  sum  in  the  denominator  of  Eq.  (2)  is  based 
on  a  test  for  equality  of  two  multinomial  distributions. 

3.  APPLICATION  TO  PHYSIOLOGICAL  DATA 

We  present  five  illustrations  of  our  approach  by 
comparing  traditional  nonlinear  measures  {D  and  K)  with 
phase-space  dissimilarity  measures  (x^  and  L).  Direct 
comparison  of  these  measures  is  difficult  due  to  their 
disparate  range,  variability,  and  physical  meaning. 
Renormalization  allows  meaningful  comparison  by 
defining  L,  as  the  value  of  each  indicator  for  the  i-th 
cutset  from  the  set,  {D,  K,  x^,  and  L}.  We  denote  the 
mean  value  of  that  indicator  over  the  base  case  by  V,  with 
a  corresponding  sample  standard  deviation,  o.  A 
renormalized  indicator  is  then  U{V)  =  |Li  -  Z  |/  c>,  which 
measures  the  number  of  standard  deviations  that  the  test 
case  deviates  from  the  baseline  mean.  Dynamical  states 
close  to  (far  from)  the  baseline  have  small  (large)  values 
of  the  renormalized  dissimilarity. 

1.  Human  electroencephalogram  data  were  acquired 
during  clinical  epilepsy  monitoring  and  analyzed  by  the 
procedure  of  Section  2.  Figure  1  shows  typical  results. 
Raw  data  in  subplot  (a)  have  very  complex,  non-periodic 
features  that  are  typical  of  brain  waves.  The  seizure  event 


occurred  at  1 10.7  minutes,  as  denoted  by  the  solid  vertical 
line  in  subplots  (d)  and  (e).  No  seizure  event  forewarning 
is  provided  by  the  correlation  dimension  in  subplot  (b),  or 
by  the  Kolmogorov  entropy  in  subplot  (c).  The  isolated 
peaks  at  42  and  58  minutes  in  subplot  (c)  are  not 
significant.  An  event  forewarning  of  27  minutes  is 
provided  by  U('^)  in  subplot  (d)  and  U(L)  in  subplot  (e), 
with  two  (or  more)  successive  occurrences  above  the 
threshold  of  5  (dashed  horizontal  line)  at  85  minutes 
(vertical  dashed  line).  Hively  and  Protopopescu  (2003) 
give  additional  details  on  the  methodology  for  this  and 
subsequent  examples. 

2.  Human  electrocardiogram  data  were  acquired 
during  ambulatory  monitoring.  Figure  2  shows  results  for 
a  ventricular  fibrillation  event  at  37  minutes.  The  raw 
data  in  subplot  (a)  show  ten  distinct  heartbeats  and  their 
associated  quasi-periodic  (nonlinear)  features.  The 
correlation  dimension  in  subplot  (b)  varies  randomly  (no 
forewarning  features)  with  a  rise  at  the  fibrillation  event. 
The  Kolmogorov  entropy  in  subplot  (c)  varies  erratically; 
the  peaks  occurring  at  8  and  24  minutes  are  not  valid 
forewarning  indications.  Event  forewarning  of  16  minutes 
(the  vertical  dashed  line)  is  provided  by  both  U('^)  in 
subplot  (d)  and  U(L)  in  subplot  (e);  forewarning 
corresponds  to  two  (or  more)  successive  occurrences 
above  the  threshold  (dashed  horizontal  line).  Similar 
results  were  obtained  for  several  additional  datasets. 

3.  Human  electrocardiogram  data  were  acquired 
during  laboratory  tests  of  fainting  (syncope),  under  the 
following  protocol:  (i)  lying  horizontal  for  10  minutes,  (ii) 
lying  in  tilted  condition  (70°  from  horizontal)  for  40 
minutes,  and  (iii)  lying  horizontal  again  for  5  minutes. 
Figure  3  (top)  shows  non-syncopal  results  with  low  values 
of  renormalized  dissimilarity  (~10)  that  increase  very 
slowly  and  erratically  over  the  tilt  period  (slope,  A~0.06- 
0.07  per  minute).  Figure  3  (bottom)  shows  syncopal 
results  for  the  same  subject  with  much  larger 
renormalized  dissimilarity  (40-70)  that  increase  much 
more  rapidly  over  the  tilt  period  (A~0.8-l  per  minute). 
The  tilt  period  in  this  second  test  was  terminated  early 
when  the  subject  fainted.  Similar  results  are  obtained  for 
a  second  subject. 

4.  Heart  wave  data  were  obtained  via  surface  chest 
electrodes  from  anesthetized  rats  subjected  to  an  induced 
sepsis  experiment.  After  55  minutes  of  normal-state 
recording,  each  test  rat  was  exposed  to  inhaled  bacterial 
endotoxin  that  induces  an  inflammatory  response  and 
eventually  sepsis.  Figure  4  shows  sample  results.  Raw 
data  in  subplot  (a)  have  distinct  heartbeats  with  additional 
quasi-periodic  (nonlinear)  features.  No  indication  of 
condition  change  is  displayed  by  either  the  correlation 
dimension  in  subplot  (b),  or  by  the  Kolmogorov  entropy 
in  subplot  (c).  The  condition  change  is  shown  clearly  by 
both  U(y^)  in  subplot  (d)  and  U(L)  in  subplot  (e),  which 
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Figure  1:  Results  for  human  electroencephalogram  channel  5  of  dataset  #PVM006,  showing  time-serial  plots  for:  (a)  2.4 
seconds  of  raw  data  collected  at  250Hz,  (b)  correlation  dimension,  D,  (c)  Kolmogorov  entropy,  K,  (d)  U{y^),  and  (e)  U{L). 
The  phase-space  dissimilarity  measures  in  subplots  (d)  and  (e)  are  for  d  =  3,  5  =  20, 1=17,  and  after  removal  of  eye  blink 
artifacts.  Each  cutest  has  N  =  22,000  points,  corresponding  to  88  seconds.  We  have  successfully  applied  this  analysis  to 
over  sixty  human  datasets 
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Figure  2:  Results  for  human  dataset  #EC8202,  showing  time-serial  plots  for:  (a)  10  seconds  of  raw  heart  wave  data 
collected  at  250  Hz,  (b)  correlation  dimension,  D,  (c)  Kolmogorov  entropy,  K,  (d)  U{y^),  and  (e)  U{L).  The  phase-space 
dissimilarity  measures  in  subplots  (d)  and  (e)  are  for  d=  5,  S  =  2i,'k  =  27,  after  removal  of  breathing  artifacts.  Each  cutest 
had  18,000  points,  corresponding  to  72  seconds. 
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Figure  3:  Results  for  human  subject  RAY  show  U(L)  and  C/(ic)  when  no  syncope  occurs  (above  the  double  line),  in 
contrast  with  U{L)  and  U{Lc)  when  syncope  does  occur  (below  the  double  line).  No  results  for  traditional  nonlinear 
measures  are  shown,  due  to  their  insensitivity  in  the  other  examples.  The  phase-space  dissimilarity  measures  are  for  d  =2,  S 
=  2,'k=  83,  after  removal  of  breathing  artifacts.  Each  cutest  has  ^=20,000  points  (80  seconds  of  data  at  a  sampling  rate  of 
250  Hz). 
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TIME  (MIN) 

Figure  4:  Results  for  dataset  #33209V,  showing  time-serial  plots  for:  (a)  2.4  seconds  of  raw  rat  heart  wave  data  collected  at 
500  Hz,  (b)  correlation  dimension,  D,  (c)  Kolmogorov  entropy  K,  (d)  U{y^),  and  (e)  U{L).  The  phase-space  dissimilarity 
measures  in  subplots  (d)  and  (e)  are  for  d  =2,  S=  2,'k  =  80,  after  removal  of  breathing  artifacts.  Each  cutest  has  A^=20,000 
points  (40  seconds  of  data). 


remains  low  for  the  first  57  minutes,  rising  abruptly  after 
the  exposure  onset,  remaining  high  for  the  next  20 
minutes,  then  decreasing  slowly  as  the  immune  response 
fought  off  the  bioagent  effects.  This  recovery  response  is 
consistent  with  other  physiological  observations  during 
the  test  (not  shown).  The  total  true  (negative)  positive 
rate  for  (un)exposed  animals  is  (6/6)  17/17. 

5.  A  surface  stethoscope  acquired  lung  sounds  data 
during  lung  experiments  on  anesthetized  pigs.  The 
baseline  state  consisted  of  normal  breathing.  Subsequent 
test  cases  were  obtained  by  injecting  a  controlled  volume 
of  air  (in  increments  of  100  milliliters  up  to  1400 
milliliters)  in  the  space  between  the  diaphragm  and  the 
lungs,  making  breathing  increasingly  more  difficult. 
Figure  5  shows  sample  pneumothorax  results.  Raw  lung 
sounds  data  in  subplot  (a)  have  very  complex  features, 
including  quasi-periodic  heartbeats  that  are 
superimposed  on  breath-cycle  undulations.  The 
correlation  dimension  in  subplot  (b)  provides  no  clear 
indication  of  condition  change.  The  Kolmogorov  entropy 
in  subplot  (c)  likewise  varies  erratically.  Condition 
change  is  indicated  by  both  U{y^)  in  subplot  (d)  and  U{L) 
in  subplot  (e);  both  rise  to  a  plateau  of  5  over  100-500 
ml,  then  increase  to  values  larger  than  20  over  500-1300 
ml  thereby  providing  robust  forewarning  of  the  animal's 
death  at  1400  ml.  Similar  results  were  obtained  for  a 
second  animal. 

4.  CONCLUSIONS 

We  have  developed  a  model-independent  method  for 
measuring  physiological  condition  change.  First,  we 
remove  confounding  artifacts  (such  as  eye -blinks  and 
breathing)  with  a  novel  zero-phase  quadratic  fdter.  The 
artifact-fdtered  data  are  converted  into  a  statistical 
distribution  function  that  describes  the  visitation 
frequency  and  sequence  of  the  dynamical  states. 
Dissimilarity  measures  between  baseline  and  test  data 
detect  the  change  by  summing  the  absolute  values  of  the 
differences  between  distribution  functions.  The 
methodology  is  quite  general  and  applies  to  a  variety  of 
physiological  data,  as  shown  by  five  examples:  (1)  brain 
waves  for  forewarning  of  epileptic  events  (Fig.  1);  (2) 
heart  waves  for  forewarning  of  ventricular  fibrillation 
(Fig.  2),  and  for  detection  of  syncope  (Fig.  3)  and  sepsis 
(Fig.  4);  and  (3)  lung  sounds  for  detection  of  breathing 
difficulty  (Fig.  5).  The  dissimilarity  measures  have  small 
values  in  the  normal  state,  followed  by  significantly 
larger  values  above  a  “normality  threshold,”  indicating 
abnormal  dynamics.  The  results  show  that  the  phase- 
space  dissimilarity  approach  is  sensitive,  robust,  and 
timely. 

Our  new  approach  also  detects  dynamical  change  in 
various  physical  processes.  Examples  include:  detecting 
balanced  and  unbalanced  centrifugal  pump  conditions 
from  motor  power;  distinguishing  different  drilling 


conditions  from  spindle  motor  current  of  a  machining 
center;  and  predicting  failure  of  a  bellows  coupling  in  a 
rotating  drive  train  from  motor  current;  discerning  the 
difference  in  micro-cantilever  vibrations  with  and 
without  mercury  on  the  sensor;  forewarning  of  failure  in 
electrical  motors;  and  forewarning  of  failure  in  motor- 
driven  components. 

We  now  have  high-fidelity  laboratory  integration  of 
the  technology  elements  into  desktop-computer  software 
that  analyzes  noisy,  archival  data  and  provides  change 
indication.  The  analysis  is  much  faster  than  real-time, 
and  can  handle  multiple  channels.  We  have  been  granted 
6  U.S.  patents  on  the  approach  (with  two  additional 
patents  pending).  Thus,  the  technology  readiness  level  is 
5.  On-going  development  includes  a  graphical  user 
interface,  more  robust  software,  and  implementation  on  a 
hand-held  device  (e.g.,  personal  digital  assistant).  We 
expect  to  complete  these  improvements  in  2004,  which 
will  allow  qualification  for  a  technology  readiness  level 
of  6. 

Success  for  such  diverse  applications  suggests  that 
the  method  can  reliably  measure  condition  change  in 
many  nonlinear  systems.  Examples  of  such  new  potential 
applications  include  combat  readiness  of  soldiers,  new 
complex  combat  systems,  and  next-generation  system- 
of-systems  (e.g.,  unit-of-action/future  combat  systems). 
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